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We investigated the structural relaxation of myosin motor domain 
from the pre-power stroke state to the near-rigor state using molec- 
ular dynamics simulation of a coarse-grained protein model. To de- 
scribe the structural change, we propose a "dual Go-model," a vari- 
ant of the Go-like model that has two reference structures. The nu- 
cleotide dissociation process is also studied by introducing a coarse- 
grained nucleotide in the simulation. We found that the myosin 
structural relaxation toward the near-rigor conformation cannot be 
completed before the nucleotide dissociation. Moreover, the re- 
laxation and the dissociation occurred cooperatively when the nu- 
cleotide was tightly bound to the myosin head. The result suggested 
that the primary role of the nucleotide is to suppress the structural 
relaxation. 

molecular dynamics simulation — dual-Go model — coarse-grained protein — 
coarse-grained nucleotide — strain sensor 

Introduction 

The mechanism of biomolecular motors is one of the major 
topics in biophysics. Among a number of such systems that 
have been found so far, the actomyosin motor is of a particular 
interest, because it is responsible for muscle contraction and 
cellular movements in eukaryotic cells. Myosin moves unidi- 
rectionally along the actin filament using chemical energy re- 
leased by ATP hydrolysis [1, 2, 3]. It is widely recognized that 
the efficiency of this energy conversion is very high compared 
with macroscopic artificial machines, in spite of the fact that 
biomolecular motors work under a noisy environment in the 
cell. In fact, the free energy released during each ATP hydrol- 
ysis is only about 20/cbT; therefore the thermal fluctuation 
should be appreciable. Although recent progress in imaging 
and of nanomanipulation has enabled the observation of single 
molecules, the movement mechanism of the actomyosin motor 
is still not understood. 

There has been a long-standing controversy between 
the tight-coupling (lever-arm) model and the loose-coupling 
model. X-ray crystallographic studies have revealed that the 
angle of the neck domain changes relative to the motor do- 
main, depending on the nucleotide state. The "lever-arm" 
model was proposed based on these observations, in which 
the structural change of myosin is tightly coupled with the 
ATP hydrolysis cycle, and directly causing a stepwise slid- 
ing motion. It was shown, however, that the sliding distance 
of the myosin along the actin filament per ATP at the mus- 
cle contraction can be much longer than the displacement 
predicted by the lever-arm-like structural change of a single 
myosin molecule [4]. Moreover, it is questionable whether a 
material as soft as proteins can accurately switch its conforma- 
tion in the same way as a macroscopic machine under thermal 



fluctuation. 

In the "loose-coupling" model, in contrast, the structural 
change does not always correspond to a step in a one-to- 
one correspondence; the motion is intrinsically stochastic and 
thermal fluctuation is an essential ingredient for its mecha- 
nism [5] . The simplest class of models that produce the loose- 
coupling mechanism is based on a thermal ratchet, in which a 
myosin molecule is treated as a Brownian particle that moves 
along a periodic and asymmetric potential under both thermal 
noise and non-thermal perturbations [6, 7]. Although ratchet 
systems can, in fact, exhibit unidirectional flow in a noisy 
environment, a high efficiency comparable to that of the acto- 
myosin system is found difficult to achieve using only a simple 
ratchet system. Even if the ratchet models do to express some 
essence of the mechanism of the biomolecular motor, they are 
too much simplified and we should say that the connection 
with the real actomyosin system is rather vague. In particu- 
lar, since the myosin is expressed as a particle, the effect of its 
conformational change is, at best, only implicitly taken into 
account. A somewhat more realistic modelling is desirable, 
which can reflect the chain conformation. 

Recently, it was revealed by single molecule experiments 
that the chemomechanical cycle of the myosin head is con- 
trolled by a load on the actomyosin crossbridge [8, 9, 10]. The 
observation suggested that the rate of ADP release from the 
myosin head depends on the force acting on myosin; namely, 
the chemical reaction rate varies with the deformation of the 
myosin. If the myosin head indeed acts as such a "strain 
sensor," this would be reminiscent of a classical model by A. 
Huxley [If]; in this model, the myosin head is supposed to 
undergo Brownian motion and change into a tightly bound 
state to the actin filament triggered by a structure dependent 
chemical reaction. 

The relationship between structure and function of pro- 
teins has long been investigated. Thus far, mainly the static 
aspects of proteins have been considered, for example, the 
classical "lock and key" model of an enzyme. Recently, the 
role of structural fluctuations, or that of a more drastic struc- 
tural change, including "partial unfolding," on protein func- 
tions has become a subject of growing interest. Although 
there have been many experimental studies to clarify the dy- 
namical processes of protein at a functional level, it is still dif- 
ficult to observe the structural changes with high resolution in 
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both space and time. Computer simulations serve as possible 
alternatives. Typical computational studies treated equilib- 
rium fluctuations near a crystal structure using the all-atom 
model [12, 13] or the elastic network model [14]; although 
these types of simulations cannot deal with large-scale struc- 
tural change, the low-frequency fluctuation modes were found 
to be consistent with the direction of motion of the structural 
change associated with the function. 

Some attempts have been made to simulate a larger struc- 
tural change beyond the elastic regime using a class of models 
called the Go-like model [15] . According to the recently devel- 
oped theory of spontaneous protein folding, the protein energy 
landscape has a funnel-like global shape toward the native 
structure. The Go-like model is certainly the simplest class of 
models that can realize a funnel-like landscape [16, 17] and has 
successfully described the folding process of small proteins. It 
is, however, not suitable for the study of a change between 
two conformations, because only the interaction between the 
pairs of residues that contact each other in the native state 
are taken into account; the conformation other than the na- 
tive state becomes too unstable as a result. Thus, a model is 
desirable in which two conformations can be embedded. Here, 
we introduce a new model bearing this property as a variant 
of the Go-like model. 

In this paper, we investigate the dynamics of the myosin 
conformational change from the pre-power stroke state to 
the near-rigor state by molecular dynamics simulations of a 
coarse-grained protein. This process is called "power stroke" 
because the angle of the lever-arm changes remarkably, and 
it is considered in the lever-arm model that this structural 
change directly causes the force generation. To describe the 
structural change, we propose the "dual Go-model." The dis- 
sociation process of the nucleotide that accompanies the con- 
formational change is also involved in the simulation by in- 
troducing a coarse-grained nucleotide. To our knowledge, the 
ligand at the binding site has not been considered explicitly in 
coarse-grained protein simulations, possibly because the pri- 
mary role of ATP is considered to be the release of chemical 
energy through hydrolysis, and the excluded volume effect of 
the molecule has not been investigated. We, however, con- 
sider that the presence or absence of nucleotides in the bind- 
ing site would profoundly affect the structural fluctuation of 
the protein. Therefore, it is important to perform simulations 
including the nucleotide. 



Results 

First, we introduce the dual Go-model. While only the native 
structure is taken as a reference structure for the potential en- 
ergy function in the standard Go- like models, the dual Go- like 
model takes two reference structures, structure 1 and struc- 
ture 2, in the effective potential. For the interaction between 
"native-contact" pairs, each potential energy function has two 
minima corresponding to two reference structures. To study 
the relaxation process from structure 2 to structure 1, the min- 
imum corresponding to structure 2 is given a slightly higher 
energy than that of structure 1 to make structure 1 more 
stable. In this study, we used a model based on one of the 
C a Go- like models [18, 19], which involves local interactions 
such as bond length, bond angle, and dihedral angle interac- 
tions as well as the native-contact interaction. For these local 



interactions, we set that each potential energy function has 
two minima of the same depth corresponding to two reference 
structures. 

As reference structures, we choose X-ray crystallography 
structures of Dictyostelium discoideum myosin II: the near- 
rigor structure without nucleotide, 1Q5G [20] for structure 
1 and the pre-power stroke structure with ADP-P; analog, 
1VOM [21] for structure 2 (Fig. 1). The initial structure 
is the pre-power stroke structure with a coarse-grained nu- 
cleotide (Fig. 2) located at the nucleotide-binding site. 

Typical time courses of the distance root mean square de- 
viation (dRMSD) are shown in Fig. 3 for (a) fc p -n = 0.6 and 
(b) 0.7, where fcp-n is the strength parameter of the protein- 
nucleotide interaction. The dRMSD from the near-rigor struc- 
ture is defined as 

dRMSD = ^Z^ZTg^T^, [ 1 ] 

in which nj = |r».j| = |rj — Tj | is the distance between C a 
carbons of the ith and jth residues in the given conformation, 
and indicates their distance in the near-rigor structure. 
At the initial conformation, the dRMSD of ~ 3.9A, and de- 
creases rapidly to ~ 3A, and stays there for a while. The 
conformation eventually relaxes into the final state, dRMSD 
~ 1.5A. This final state is actually the near-rigor state, judg- 
ing from its average structure; in fact, the dRMSD of the 
average structure is ~ lA. In short, the myosin motor domain 
in the pre-power stroke conformation relaxes at first into the 
intermediate state and then to the near-rigor conformation. 

We introduced an index to characterize the state of the 
nucleotide binding, Q nuc i(r); we count how many of the nu- 
cleotide contacts that is formed in structure 2 (the pre-power 
stroke) remain in a given conformation, Y. Then, Q nuc i(r) is 
defined as this number divided by the number of nucleotide 
contacts in structure 2. Q nuc i ~ 1 when a nucleotide is bound, 
and Q nuc i = if the nucleotide-binding site is empty. The 
time courses of Q nuc i are also shown in Fig. 3. The structural 
relaxation occurs after or, at the earliest, at the same time as 
the nucleotide dissociation. Furthermore, the relaxation tends 
to synchronize with the dissociation as fcp-n increases. To clar- 
ify the fcp-n dependency of the synchronization, we plotted the 
histograms of and of At — — r r from 200 independent 
runs for each value of fcp-n, where is the number of steps 
taken before the nucleotide dissociates, r r is the number of 
steps taken before the conformation relaxes to the near-rigor 
state, and At is the delay in the relaxation after the dissoci- 
ation takes place (Fig. 4). 

For small fcp-n, both the histograms of both and of At 
show exponential decay; thus, the nucleotide dissociation and 
the relaxation of the myosin conformation are considered to be 
decoupled. For larger fcp-n, on the other hand, the histogram 
of cannot be fitted to an exponential decay. In addition, 
the average of is shifted to the right and the delays, At be- 
come shorter; in other words, the nucleotide is unbound later 
and the conformational relaxation tends to occur immediately 
after the nucleotide dissociation. For fc p - n = 0.7, dissociation 
and relaxation occur nearly simultaneously in over 60% of 200 
trajectories. Note that apparent "At < cases" are caused 
simply from the numerical ambiguity of tj and T r and actually 
correspond to coincidental dissociation-relaxation. 

The largest difference between the intermediate state and 
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the initial (pre-power stroke) conformation is the position of 
the converter domain with respect to the other subdomains; 
the relative positions among other subdomains (for exam- 
ple, the N terminal and the 50-kDa subdomain) are similar 
to those in the pre-power stroke conformation. The average 
dRMSD of the intermediate state varies slightly with the pa- 
rameter kp-n (Fig. 5), reflecting little difference in the position 
of the converter relative to the other subdomains. 

Some of the native contacts of structure f are not formed 
until the conformation finally relaxes to the near-rigor state. 
Figure 6 shows the residues included in these contacts. They 
are concentrated at the boundary between the N terminal 
and the 50-kDa subdomain, that is, the region around the 
nucleotide-binding site. Thus, the final relaxation process con- 
sists of a rearrangement of the N terminal subdomain against 
the other part. 

As already mentioned, the myosin motor domain relaxes to 
the near-rigor conformation only after the dissociation of the 
nucleotide and not before. Thus, it seems that the nucleotide 
must be unbound for the final relaxation to occur. This ob- 
servation leads to a speculation that the nucleotide blocks the 
deformation of myosin around the nucleotide-binding site by 
its volume. To investigate the case of when the nucleotide 
cannot dissociate, we attempted a nucleotide-free and con- 
strained simulation. In this simulation, instead of treating the 
nucleotide molecule explicitly, we connected the residues that 
would contact with the nucleotide by virtual bonds, to force 
the nucleotide-binding site to keep the pre-power stroke form. 
Figure 5 shows the time courses obtained by the simulations. 
We find that the conformation remains at the intermediate 
state and that the relaxation toward the near-rigor state is 
not completed. 

Discussion 

The relaxation simulation using coarse-grained myosin and 
the nucleotide have shown that the myosin motor domain does 
not relax to the near-rigor conformation before the nucleotide 
dissociates. Ishijima el al. [22] showed by simultaneous ob- 
servation of ADP release and mechanical events that force is 
generated at the same time as or several hundreds of millisec- 
onds after the dissociation of ADP. Our results are consistent 
with their experimental findings if force generation is preceded 
by structural relaxation. Moreover, the results from the sim- 
ulations in which the conformation of the nucleotide-binding 
site is constrained also indicate that the relaxation is indeed 
prevented if the nucleotide cannot dissociate. 

Based on these observations, we now suggest that the pri- 
mary role of the nucleotide in the "power stroke" process is 
to suppress relaxation through blocking deformation around 
the nucleotide-binding site by its volume. In this scenario, hy- 
drolysis is required to alter the affinity of the nucleotide to the 
binding site. In particular, our simulations have shown that 
the structural relaxation is synchronous with nucleotide dis- 
sociation when the nucleotide is tightly bound to the myosin 
head. In other words, the nucleotide dissociates cooperatively 
with the motion of the subdomain. This strong coupling of 
deformation and dissociation seems to be relevant to the func- 
tion of the "strain sensor," in which nucleotide dissociation is 
controlled by the strain induced by an external force. The cor- 
relation depends on the binding strength, fcp-n; the relaxation 
is only loosely coupled with the dissociation for weak bind- 
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ing conditions. The origin of a large kinetic diversity among 
myosins [23, 24] may be attributed to this binding strength 
dependence. 

The intermediate state observed in the relaxation process 
should also be discussed. Although several intermediate states 
have been revealed by kinetic experiments[23], their struc- 
tural aspects, except for ADP-P; state, are little known. Shih 
et al. [25] reported, from their FRET study, that there are 
two "pre-power stroke" conformations; while one conforma- 
tion corresponds to the crystal structure of the complex with 
the ADP-Pj analog, the other conformation has not yet been 
observed using crystallography. We found that the average 
structure of the intermediate state observed in the present 
study is consistent with the latter conformation. 

Our dual Go model is effective in studying myosin con- 
formational changes. Recently, a similar approach was pro- 
posed for a Monte Carlo simulation of a lattice protein, in 
which a double-square-well potential for native-contact pairs 
was introduced [26]. This type of model, in which two confor- 
mations are embedded in an energy potential function, will 
be useful in understanding the dynamics of protein conforma- 
tional change. 

Models and Methods 

Our dual Go-model is a variant of the C a Go-like model 
[18, 19]. A protein chain is formed, consisting of spherical 
beads that represent C a atoms of amino acid residues con- 
nected by virtual bonds. In conventional Go-like models, only 
amino acid pairs that contact in the native conformation are 
assigned an effective energy. In the dual Go-model, on the 
other hand, the effective energy function takes two reference 
structures, structure 1 and structure 2. A nucleotide molecule 
is also expressed as a chain of connected beads. 

The total energy of the system, <7 tot consists of three 
terms; 

U tot = Up + Un + Up-n, [2] 

where Up is the intraprotein interaction, Un is the intranu- 
cleotide interaction, and Up-n is the interaction between pro- 
tein and nucleotide. 

The effective protein energy Up at a conformation T is 
given as, 

[/ P (r,r (1) ,r (2) ) = u b + u e + u+ + u nc + u nnc . [3] 

where T^ 1 - 1 and stand for the conformations of the two 
reference structures. The terms in Eq. 3 are defined as fol- 
lows: 

,7 z = ^min{\/ z (1) ( 2l ),V z (2) (^)}, [4] 

i 

native 

U nC = C °ff min{yW(r«), Ci^fo)}, [5] 

j<i-3 

non-native 
contact 

C/ mC = £ VS nC (^), [6] 

j<i-3 

where z stands for b, 6 or (j>, and (1) and (2) again indicate 
the reference conformations. The vector Vij = r; — Yj is the 
distance between the ith and jth C a , where r; is the position 
of the ith Co,. 6j = |b»| = |r»j + i| is the virtual bond length 
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between two adjacent C Q . 8i is the angle between two adja- 
cent virtual bonds, where cosf?i = bi_i ■ b;/6j_i&i, and <pi is 
the jth dihedral angle around hi. The first three terms of Eq. 
3 provide local interactions, while the last two terms are in- 
teractions between non-local pairs that are distant along the 
chain. 



(a) T/ («) 



For potential functions, V 2 
same functions as Clementi el al. [18]: 



V^\b i )=k h {b i -br ) ) 



(a),2 



[7] 



g(«)\2 



V, 



(«)/ 



1 — cos( 



(0i-^ a) )) 
±(l-coe3(&-^ a) )) 



Vnc^Hj) = fcnc 



,(«) 



^j mC ( r y) = fe nnc ( — 



[101 



[11] 



where the superscript a is 1 or 2 and represents the appro- 
priate reference structure. Parameters with superscript 1 or 
2 are constants taken from the corresponding values in struc- 
ture 1 or 2, respectively. For the local interaction terms (bond 
length, bond angle, and dihedral angle), the potential energy 
for each set of beads takes the smaller of and V' 2 '. For 
example, the length of the ith bond is b' 1 ' in structure 1 and 
is in structure 2; therefore, the potential energy for this 
bond is fc b min{(& l - bf ) ) 2 , (b t - bf ] ) 2 }. We specify(define) 
that the jth and jth amino acids are in the "native contact 
pair" of structure 1 (or structure 2) if one of the non-hydrogen 
atoms in the jth amino acid is within 6.5 A of one of the non- 
hydrogen atoms in the jth amino acid at the structure 1 (or 
structure 2). The interaction potential for each native-contact 



pair, ij, takes the smaller of V^c ( r »i) an d Ci2Vri c (r^). C12 is 
the ratio of the potential depth of structure 2 to that of struc- 
ture 1 (Fig. 7). Here, since we intend to perform simulations of 
the structural change from structure 2 to structure 1, that is, 
we want to structure 1 as the final stable structure, we assign 
a value smaller than unity to C12 (C12 = 0.8). If a residue pair 
ij is a native-contact pair in structure 1 but not in structure 2 
(or vice versa), the interaction potential between the ith and 
jth residues is a Lennard- Jones potential, Vnc^ or2 ^ (r 4 j), with 
a single minimum. Other relevant parameters are — 100.0, 
k = 20.0, k<f, = 1.0, k nc = fcnnc = 0.25, and C = 4.0. The 



cut-off length for calculating ' is taken to be 2r\ . 

To understand the mechanism of the actomyosin motor, 
it is desirable to study the conformational change to the rigor 
state. However, currently no X-ray crystal structure is cur- 
rently available for a true rigor complex with actin; there- 
fore, we studied the structural relaxation from the pre-power 
stroke state to the near-rigor state. Structures 1 and 2 thus 
correspond to the near-rigor and pre-power stroke structures, 
respectively. We used 1Q5G [20], which is the nucleotide 
free structure of Dictyostelium discoideum myosin II, as the 
near-rigor state. Although a few nucleotide- free structures of 



» 



myosin II have been determined, only 1Q5G is regarded to be 
the near-rigor state, because both switch I and switch II are 
in the open position. We chose 1VOM [21] for the pre-power 
stroke structure. It includes an ADP-Pj analog (ADP-VCU) in 
the nucleotide binding site. While 1Q5G consists of residues 2- 
765 without a gap region, 1VOM includes only residues 2-747 
and has gap regions where the structure has not been deter- 
mined by X-ray crystallography. Therefore we used the struc- 
tures of residues 2-747 for simulations. Potential functions of 
the near-rigor conformation are assigned for local interactions 
in the gap region. 

The nucleotide molecule (ADP-V0 4 in 1VOM) is also in- 
cluded in the simulation as a coarse-grained chain (Fig. 2). 
The coarse-grained ADP-V04 is represented as a short linear 
chain of five beads, corresponding to a purine base, a sugar 
(ribose), two phosphates, and VO4 (a phosphate analog). The 
intranucleotide interaction is defined as 



Un 



r k +i - r k \ 



r (2) 

l k + l 



(2) h 2 



+E fc b(i 



ffc+2 — rfc| 



r (2) 



ri 2) |) 2 . [12] 



The interaction potential between the protein and the nu- 
cleotide is similar to that between the non-local residues in the 
protein. Here, we assume that only structure 2 includes the 
nucleotide: therefore, the potential function has only a single 
well (the standard Lennard- Jones potential). We specify that 
the ith residue of the protein and the fcth bead in the nu- 
cleotide chain should be in "native-contact" in the pre-power 
stroke conformation when one of the non-hydrogen atoms in 
the kth bead (base or sugar of P;) is within 4.5 A of one of the 
non-hydrogen atoms in the jth amino acid. The residues that 
form native-contacts with the nucleotide are called nucleotide- 
contact residues: 



native 
c onta ct 

Up-n(r ik ) = ^ kp-n 

i,k 

non-native 
contact 

+ E ' 

3<i-3 



.(2) 



Tik 



-6 



r (2)~ 
' ik 

Tik 



[13] 



where i stands for the ith residue and k stands for the feth 
bead in the nucleotide chain. 

The dynamics of the proteins are simulated using the 
Langevin equation at a constant temperature T, 



[141 



where v is the velocity of the ith bead and a dot represents the 
derivative with respect to time t (thus, v» = r^), and Fj and 
£» are systematic and random forces acting on the iih bead, 
respectively. The systematic force F; is derived from the effec- 
tive potential energy U and can be defined as F; = —dU/dri. 
£i is a Gaussian white random force, which satisfies = 
and = 2fT5ij5(t — t')l, where the bracket denotes 

the ensemble average and 1 is a 3 x 3 unit matrix. We used 
an algorithm by Honeycutt and Thirumalai [27] for a numer- 
ical integration of the Langevin equation, We used 7 = 0.25, 
irii=1.0, and the finite time step At = 0.02. 

For a given protein conformation, F, we note that the na- 
tive contact of structure 1 (or 2) between i and j is formed 
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if the C a distance Vij = \iij\ satisfies 0.8r^ < rij < 1.2r^ 

(0.8rg> <ry < 1.2r<f). 

Simulations were started from the pre-power stroke struc- 
ture. The initial positions of residues in the gap regions of 
1VOM were set randomly under the condition that the bond 
length was 3.8 A. The initial velocities of each bead was given 
to satisfy the Maxwell distribution. The temperature was set 
lower than the folding temperature for structure 1. 

We also ran a nucleotide-free and constrained simulation, 
in which the nucleotide was not explicitly included but the 
relative positions of the nucleotide-contact residues were con- 
strained by virtual bonds in all-to-all correspondence to keep 
the pre-power stroke form. The natural length of the virtual 



(2) 

bonds i-j, r\ • , is the C a distance between the ith and jth 
residues in the pre-power stroke conformation (structure 2). 
The total effective potential energy was Up + Ucon , and 

nucl.- 
c onta ct 

C/con= J2 k C on( rij -4f) 2 , [15] 

3<i 

where fc C on = 1.0. 
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Fig. 1. We chose (a) near-rigor structure, 1Q5G, and (b) pre-power stroke structure, 1V0M for structure 1 and 2, respectively. 1V0M contains the ADP-Pj analog, 
ADP-V04 (yellow). Also shown are the N-terminal (green), 50 kDa subdomain (red) and the converter (cyan) included in C-terminal subdomain (blue) that is connected to 
the lever arm. 




Fig. 2. (a) ATP and (b) coarse-grained ATP. 
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Fig. 4. Histogram of (a) the number of steps before dissociation and (b) the delay of the relaxation after the dissociation from 200 independent runs for each kp-n- 
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Fig. 6. Residues included in the contacts that are formed at the final relaxation to the near-rigor. 
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